## Script to run model and plots of rebel votes by Wordscores extremism
## Replicates Table 3 and Figure 2 in paper
## JBS: final revision revised 3 August 2017

rm(list=ls())

library(lme4)
library(stargazer)
library(mvtnorm)
library(arm)

##### SET YOUR PATH HERE ##############
##########################################

path <- "~/Dropbox (Personal)/Rebel Summaries/APSR_SKLLO_Repfiles/"

setwd(path)

# Load Wordscores Positions for Labour
load("CalcData/positions_jsnew.Rda")

# Load rebellion data
load("CalcData/rebeldta0510.Rda")
load("CalcData/rebeldta9297.Rda")

## Start with Labour
# Keep only pre-Blair scores

ws.pos <- results1[results1$firstperiod==1,]

# Keep only Labour for now
reblab92 <- all.dta.9297[all.dta.9297$Party=="Lab",]

fname <- as.character(reblab92$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(reblab92$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

reblab92$fullname <- paste(fname, lname, sep=" ")

## Fix names that differ

# Changed Alan Williams from Swansea because it's A. Williams from Carmarthen who should match
reblab92$fullname[reblab92$MPName.x=="Alan John Williams"] <- "Alan John Williams"

ws.pos$names[ws.pos$names=="Tony Benn"] <- "Anthony Benn"
ws.pos$names[ws.pos$names=="Eddie O'Hara"] <- "Edward O'Hara"
ws.pos$names[ws.pos$names=="Frank Cook"] <- "Francis Cook"
ws.pos$names[ws.pos$names=="Jeff Rooker"] <- "Jeffrey Rooker"
ws.pos$names[ws.pos$names=="Ronnie Campbell"] <- "Ronald Campbell"
ws.pos$names[ws.pos$names=="Terry Rooney"] <- "Terrence Rooney"
ws.pos$names[ws.pos$names=="Tom Clarke"] <- "Thomas Clarke"
ws.pos$names[ws.pos$names=="Tony Wright"] <- "Anthony Wright"
ws.pos$names[ws.pos$names=="Win Griffiths"] <- "Winston Griffiths"
ws.pos$names[ws.pos$names=="Jim Cousins"] <- "James Cousins"
ws.pos$names[ws.pos$names=="Jim Cunningham"] <- "James Cunningham"
ws.pos$names[ws.pos$names=="Jimmy Wray"] <- "James Wray"
ws.pos$names[ws.pos$names=="Jo Richardson"] <- "Josephine Richardson"
ws.pos$names[ws.pos$names=="Bill Michie"] <- "William Michie"
ws.pos$names[ws.pos$names=="Eddie Loyden"] <- "Edward Loyden"
ws.pos$names[ws.pos$names=="Harry Barnes"] <- "Harold Barnes"
ws.pos$names[ws.pos$names=="Ken Livingstone"] <- "Kenneth Livingstone"
ws.pos$names[ws.pos$names=="Llin Golding"] <- "Llinos Golding"
ws.pos$names[ws.pos$names=="Mo Mowlam"] <- "Marjorie Mowlam"
ws.pos$names[ws.pos$names=="Ron Davies"] <- "Ronald Davies"
ws.pos$names[ws.pos$names=="Ron Leighton"] <- "Ronald Leighton"
ws.pos$names[ws.pos$names=="Sam Galbraith"] <- "Samuel Galbraith"
ws.pos$names[ws.pos$names=="Ted Rowlands"] <- "Edward Rowlands"

mergedta.92 <- merge(reblab92 , ws.pos, by.x = "fullname", by.y = "names")

mergedta.92$Score <- -1*mergedta.92$Score

### Later period

ws.pos <- results1[results1$firstperiod==0,]

# Keep only Labour for now
reblab05 <- all.dta.0510[all.dta.0510$Party=="Lab",]

# Load rebel data and create name var for merging

fname <- as.character(reblab05$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(reblab05$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

reblab05$fullname <- paste(fname, lname, sep=" ")

## Fix names
## Changed Alan Williams. Only one Alan Williams in 2nd period

ws.pos$names[ws.pos$names=="Alan WilliamsSwansea"] <- "Alan Williams"
ws.pos$names[ws.pos$names=="Huw IrrancaDavies"] <- "Huw Irranca-Davies"
ws.pos$names[ws.pos$names=="Robert MarshallAndrews"] <- "Robert Marshall-Andrews"

mergedta05 <- merge(reblab05 , ws.pos, by.x = "fullname", by.y = "names")

mergedta05$Score <- -1*mergedta05$Score

### Models for Labour 92 and 05
## Models 9 and 11 in Table 3
mod1 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta.92)
summary(mod1)

mod2 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta05)
summary(mod2)


####################
## Conservatives
####################

## remove objects used Labour estimation
rm(list=c("fname","lname","mergedta.92","mergedta05","reblab05","reblab92","results1","ws.pos" ))

## Load Wordscores positions for Conservatives
load("CalcData/positions_tories.Rda")

# Keep only first period scores

ws.pos <- results1[results1$firstperiod==1,]

# Keep only Cons 
rebcon92 <- all.dta.9297[all.dta.9297$Party=="Con",]

fname <- as.character(rebcon92$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(rebcon92$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

rebcon92$fullname <- paste(fname, lname, sep=" ")

## Fix names that differ
ws.pos$names[ws.pos$names=="Bill Cash"] <- "William Cash"
ws.pos$names[ws.pos$names=="Bob Dunn"] <- "Robert Dunn"
ws.pos$names[ws.pos$names=="Bob Spink"] <- "Robert Spink"
ws.pos$names[ws.pos$names=="Charles GoodsonWickes"] <- "Charles Goodson-Wickes"
ws.pos$names[ws.pos$names=="Elaine KellettBowman"] <- "Dame Kellett-Bowman"
ws.pos$names[ws.pos$names=="Jill Knight"] <- "Dame Knight"
ws.pos$names[ws.pos$names=="Iain DuncanSmith"] <- "Iain Duncan Smith"
ws.pos$names[ws.pos$names=="James DouglasHamilton"] <- "Lord Douglas-Hamilton"
ws.pos$names[ws.pos$names=="Michael McNairWilson"] <- "Patrick McNair-Wilson"
ws.pos$names[ws.pos$names=="Tom King"] <- "Thomas King"
ws.pos$names[ws.pos$names=="Tony Baldry"] <- "Anthony Baldry"
ws.pos$names[ws.pos$names=="Tony Marlow"] <- "Antony Marlow"
ws.pos$names[ws.pos$names=="Tim Boswell"] <- "Timothy Boswell"
ws.pos$names[ws.pos$names=="Tim Devlin"] <- "Timothy Devlin"
ws.pos$names[ws.pos$names=="Tim Smith"] <- "Timothy Smith"

mergedta.92 <- merge(rebcon92 , ws.pos, by.x = "fullname", by.y = "names")

## Later period

ws.pos <- results1[results1$firstperiod==0,]

# Keep only Cons

rebcon05 <- all.dta.0510[all.dta.0510$Party=="Con",]

# Load rebel data and create name var for merging

fname <- as.character(rebcon05$MPName)

fname <- sapply(strsplit(fname, split=" "), "[" , 1)

lname <- as.character(rebcon05$Name)

lname <- sapply(strsplit(lname, split=","), "[" , 1)

rebcon05$fullname <- paste(fname, lname, sep=" ")


## Fix names
## Changed Alan Williams. Only one Alan Williams in 2nd period

ws.pos$names[ws.pos$names=="David HeathcoatAmory"] <- "David Heathcoat-Amory"
ws.pos$names[ws.pos$names=="Geoffrey CliftonBrown"] <- "Geoffrey Clifton-Brown"
ws.pos$names[ws.pos$names=="Ian LiddellGrainger"] <- "Ian Liddell-Grainger"
ws.pos$names[ws.pos$names=="Ian LiddellGrainger"] <- "Ian Liddell-Grainger"

mergedta05 <- merge(rebcon05 , ws.pos, by.x = "fullname", by.y = "names",all.y=T)


### Models for Cons 92 and 05
## Models 10 and 12 Table 3
mod3 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta.92)
summary(mod3)

mod4 <- glmer(matrix(c(tot.reb , tot.withparty), ncol=2) ~ gov*Score+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), control=glmerControl(optimizer="bobyqa"), data= mergedta05)
summary(mod4)

stargazer(mod1,mod3,mod2,mod4,type="latex",title="Effect of Government Status on Rebellion: Fractional Logit Models with Random Effects (Extremism Scores)", align=T, column.labels=c("Labour","Conservative","Labour","Conservative"),dep.var.labels=c("1992--2001","2005--2015"),style="ajps")

#### Create Plots of Simulations
#### Simulate backbenchers with mean victory margin (majority)

labdta9297 <- all.dta.9297[all.dta.9297$Party=="Lab",]
condta9297 <- all.dta.9297[all.dta.9297$Party=="Con",]
labdta0510 <- all.dta.0510[all.dta.0510$Party=="Lab",]
condta0510 <- all.dta.0510[all.dta.0510$Party=="Con",]

# Set number of simulations
sims<-1000

# Take random draws
dist<-rmvnorm(sims, fixef(mod1), as.matrix(vcov(mod1))) # Labour 92
dist2<-rmvnorm(sims, fixef(mod2), as.matrix(vcov(mod2))) # Labour 05
dist3<-rmvnorm(sims, fixef(mod3), as.matrix(vcov(mod3))) # Con 92 
dist4<-rmvnorm(sims, fixef(mod4), as.matrix(vcov(mod4))) # Con 05

# Z is indicator of in gov or not
z<-c(0,1)

calc<-length(z)

lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
lo5<-numeric(calc)
hi5<-numeric(calc)
pe5<-numeric(calc)
lo6<-numeric(calc)
hi6<-numeric(calc)
pe6<-numeric(calc)
lo7<-numeric(calc)
hi7<-numeric(calc)
pe7<-numeric(calc)
lo8<-numeric(calc)
hi8<-numeric(calc)
pe8<-numeric(calc)

for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  pp5<-numeric(sims)
  pp6<-numeric(sims)
  pp7<-numeric(sims)
  pp8<-numeric(sims)
  
  for(j in 1:sims){
    #Labour (92-97); Moderate
    pp[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,7]*0*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T) + dist[j,6]*mean(labdta9297$start_year,na.rm=T)
    #Labour (92-97); extreme
    pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,7]*1*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T) + dist[j,6]*mean(labdta9297$start_year,na.rm=T)
    #Conservative (92-97); Moderate
    pp3[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*0+dist3[j,7]*0*z[i]+dist3[j,4]*mean(condta9297$Majority, na.rm=T) + dist3[j,6]*mean(condta9297$start_year,na.rm=T)
    #Conservative (92-97); Extreme
    pp4[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*1+dist3[j,7]*1*z[i]+dist3[j,4]*mean(condta9297$Majority, na.rm=T) + dist3[j,6]*mean(condta9297$start_year,na.rm=T)
    #Labour (05-10); Moderate
    pp5[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*0+dist2[j,7]*0*z[i]+dist2[j,4]*mean(labdta0510$Majority, na.rm=T) + dist2[j,6]*mean(labdta0510$start_year,na.rm=T)
    #Labour (05-10); extreme
    pp6[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*1+dist2[j,7]*1*z[i]+dist2[j,4]*mean(labdta0510$Majority, na.rm=T) + dist2[j,6]*mean(labdta0510$start_year,na.rm=T)
    #Conservative (05-10); Moderate
    pp7[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*0+dist4[j,7]*0*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T) + dist4[j,6]*mean(condta0510$start_year,na.rm=T)
    #Conservative (05-10); Extreme
    pp8[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*1+dist4[j,7]*1*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T) + dist4[j,6]*mean(condta0510$start_year,na.rm=T)
  }
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  pe5[i]<-invlogit(mean(pp5))
  lo5[i]<-invlogit(quantile(pp5, 0.025))
  hi5[i]<-invlogit(quantile(pp5, 0.975))
  
  pe6[i]<-invlogit(mean(pp6))
  lo6[i]<-invlogit(quantile(pp6, 0.025))
  hi6[i]<-invlogit(quantile(pp6, 0.975))
  
  pe7[i]<-invlogit(mean(pp7))
  lo7[i]<-invlogit(quantile(pp7, 0.025))
  hi7[i]<-invlogit(quantile(pp7, 0.975))
  
  pe8[i]<-invlogit(mean(pp8))
  lo8[i]<-invlogit(quantile(pp8, 0.025))
  hi8[i]<-invlogit(quantile(pp8, 0.975))
  
  print(i)
}


#92-2001
pdf("FinalPlots/Reb_by_Ext_92_2001.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe1, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.17), ylab="Probability of Rebellion", xaxt="n", col="grey80", xlab="")
segments(0.15, lo1[1], 0.15, hi1[1], lwd=2, col="grey80")
segments(0.65, lo1[2], 0.65, hi1[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe2, pch=16, cex=1.5, col="grey80")
segments(0.20, lo2[1], 0.20, hi2[1], lwd=2, col="grey80")
segments(0.70, lo2[2], 0.70, hi2[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe3, pch=15, cex=1.5, col="grey20")
segments(0.25, lo3[1], 0.25, hi3[1], lwd=2, col="grey20")
segments(0.75, lo3[2], 0.75, hi3[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe4, pch=16, cex=1.75, col="grey20")
segments(0.30, lo4[1], 0.30, hi4[1], lwd=2, col="grey20")
segments(0.80, lo4[2], 0.80, hi4[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Moderates", "Extremists", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()

#2005-2010
pdf("FinalPlots/Reb_by_Ext_05_2015.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe5, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.17), ylab="Probability of Rebellion", xaxt="n", col="grey80", xlab="")
segments(0.15, lo5[1], 0.15, hi5[1], lwd=2, col="grey80")
segments(0.65, lo5[2], 0.65, hi5[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe6, pch=16, cex=1.5, col="grey80")
segments(0.20, lo6[1], 0.20, hi6[1], lwd=2, col="grey80")
segments(0.70, lo6[2], 0.70, hi6[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe7, pch=15, cex=1.5, col="grey20")
segments(0.25, lo7[1], 0.25, hi7[1], lwd=2, col="grey20")
segments(0.75, lo7[2], 0.75, hi7[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe8, pch=16, cex=1.75, col="grey20")
segments(0.30, lo8[1], 0.30, hi8[1], lwd=2, col="grey20")
segments(0.80, lo8[2], 0.80, hi8[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Moderates", "Extremists", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()




